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V . 
f^ We propose a simplification of the Optimized Effective Potential (OEP) method 

applied to the Self Interaction Correction (SIC) on the Local Density Approximation 

(LDA) in Density Functional Theory (DFT). The new scheme fulfills crucial formal 

key properties. It turns out to be simple and accurate. We apply the new method to 

a schematic model for a dimer molecule and to the C atom. We discuss observables 

which are especially sensitive to details of the SIC. 
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5 ■ 1 Introduction 

O 



Density-functional theory (DFT) is a major tool for the description of elec- 
tronic systems as, e.g., atoms, molecules, solids, or chemical reactions [I]. It 
simplifies the involved many-electron problem and thus allows to compute 
rather complex systems, for a general overview sec [2.3J. Most applications 
employ energy-density functionals based on the Local Density Approximation 
(LDA), see e.g. [3], or its extension to the generalized gradient approximation 
(GGA) [5]. In spite of their success, these approaches have still deficiencies. 
In particular, the self-interaction error spoils single-particle properties as, e.g., 
the Ionization Potential (IP) or the band gap in solids PIT] . Another critical 
detail where LDA and GGA fail is the polarizability in chain molecules [8119] . 
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A way out of the self- interaction problem is to use exact exchange. That, how- 
ever, is an order of magnitude more expensive and it causes problems with a 
reliable description of the remaining correlation effects [TU] . An intuitive and 
efficient solution, still close to the spirit of DFT, is to augment LDA by a 
Self- Interaction Correction (SIC) [TTfT2] . That, on the other hand, produces 
a state-dependent mean-field Hamiltonian which requires extra measures to 
maintain orthogonality of the single particle basis [T2|13|I14|[T5] . The desirable 
goal remains to have a common local mean-field potential V(r). That goal is 
reached through the optimized effective potential (OEP) method, see [TSlfTT] 
for early proposals, and [TU] for an extensive recent review. There exist OEP 
approaches based on exact exchange but these are still rather involved, as 
exchange is. In that case, OEP serves mainly to produce a well defined single- 
particle spectrum in occupied space as well as in unoccupied space [T8][T9] . 
But OEP is also a useful means to overcome the difficulties of SIC, see e.g. 
[2TJ|21f22| . Indeed it is found that OEP manages to maintain crucial features 
of the underlying SIC (or exact exchange) as, e.g., the localization of states or 
the derivative discontinuity [23|24|25j . The involved SIC-OEP equations are 
thus often simplified. Popular is, e.g., the Krieger-Li-Iafrate (KLI) approach 
[2"U]r2"T] and, in a further step of simplification, the Slater approximation |16j . 
However, KLI and Slater approximation can easily miss crucial features of 
SIC as, e.g., the localization of states and the performance with respect to 
polarizability [8f9] . 



It is the aim of this paper to propose a scheme for SIC-OEP which allows to 
deal with the Slater approximation to OEP while maintaining key features of 
the full SIC scheme. The new move is to allow for a double set of occupied 
single-particle states. The two sets are connected by a unitary transformation 
such that both sets build the same total density. The first of the sets consists 
of the solutions of the SIC-OEP equations and is thus diagonal in energy. The 
second set is used to compute the total SIC energy. The distinction between 
the two sets allows the second set to develop spatially localized states which 
is energetically advantageous for the SIC energy functional. The localization 
of the second set validates the Slater approximation to OEP. We call the new 
scheme "generalized Slater approximation" . We present it after briefly review- 
ing SIC-OEP, SIC as such, and a variant of SIC which already employs the 
double set of single-particle states. The method is applied to a schematic model 
for a dimer molecule and to the C atom. We will compare exact Hartree-Fock 
with all relevant approximations (LDA, SIC, Slater and generalized Slater). In 
all the following, we employ, in fact, the local spin-density approximation but 
denote it with the simpler term LDA (rather than the sometimes used more 
complete abbreviation LSDA). 



2 SIC without approximations 



The starting point for all following considerations is the SIC energy functional 
for electrons, which reads 



N 

EsiC = Ekin + E ion + -ExDA \fi] — X! -^LDA [pfi] ( 1 ) 

(5=1 

where £"lda[p] is a standard energy-density functional in LDA complementing 
the kinetic energy E^ n and the interaction energy with the ionic background 
E- lon . The last term therein is the SIC term. The densities are defined from 
the single-particle states if)p, i.e. p = J2p Pp with p@ = l^l 2 - All summations 
run over occupied states only. 

Two features of the SIC energy (CQ) are to be noted. First, it is not invariant 
under unitary transformations amongst the occupied states because the SIC 
term is sensitive to each single-particle state separately. Second, the SIC in- 
duces a tendency to localized single-particle states. Given a total density p, the 
distribution over several pp produces maximal Coulomb energy if these pp are 
localized. Such a localization is found in many studies of SIC [26] or SIC-OEP 
[ID] . This feature will play a crucial role in the following considerations. 

The mean-field equations are derived from the total energy (Tj[|) by variation 
with respect to the single-particle states ^* . In the general case, the mean-field 
Hamiltonian becomes 



h a = ti LD ^-U a (v) , (2a) 

U a (r) = S -^^ = U LDA [\U] , (2b) 

^(LDA) = £_ + ^DAjp] (2c) 

2m op(r) 

where the first term is the standard LDA mean field and the second term, 
state-dependent, stems from the SIC term in the energy (fll). 



2.1 SIC with single set 



The SIC mean- field Hamiltonian ([2]) is state-dependent. A direct solution of 
the Schrodinger equation with a state-dependent Hamiltonian can lead to non- 
orthogonal single-particle states. This, however, violates basic requirements 
of DFT. One has to enforce orthonormality by a constraint, thus minimizing 



Esic — J2af3 'W^alV'/?) where \ a p is a matrix of Lagrangian multipliers, which 
is non-diagonal in general. By straightforward variational techniques, we then 
obtain the stationary SIC equations : 



^sic#a) = ]Cl^i)V* > (3a) 

^p a = {ipp\hs\c\^o) , (3b) 

Q={il)p\Up-U a \il) a ), (3c) 

h$ic = h- L v A -^U a \'4) a ){'4) a \ , (3d) 

a 

where we have packed the state-dependence into projector notation in the 
SIC mean field (13dl) . This yields a formally compact notation of h$ic- Note 
that Eq. fl3a| is not an eigenvalue equation because the matrix of Lagrange 
multipliers may possibly become non-diagonal. This equation as such thus does 
not provide a valuable starting basis for developing an OEP scheme, which 
requires implicitly an eigenvalue equation [T7|. A key part in the SIC equations 
([3]) is the symmetry condition (l3c]) . It emerges as follows : The Lagrangian 
matrix is real symmetric, X a p = Xp a , because the overlap matrix (4> a \ip/3) 
has that feature. Thus the hermitian conjugate of (I3"a| becomes (il>p\h' SIC = 
J2 a -Va ("?/>« | • We take the difference. The \p a cancel out and we remain with 
= C0,a|^sic ~~ ^sicI^q) which then yields the symmetry condition (I3cl) . It 
is a highly non-linear equation requiring an involved solution strategy. But 
formally, the symmetry condition simply emerges from minimizing the total 
energy in a reduced space of orthonormalized single-particle orbitals. This 
minimization principle also implies that there always exists a solution to the 
symmetry condition, in spite of its apparent complexity. 



2.2 SIC with double basis set 



The SIC equations ([3]) do not yield immediately single-particle energies and 
associated "energy- diagonal" states. It is plausible that the matrix \p a con- 
tains the relevant information. The idea is thus that diagonalization of X/3 a 
yields the desired single-particle energies together with energy-diagonal states 
(fii. These energy- diagonal states tend naturally to larger spatial width and 
will not remain as localized as the SIC basis states ip a . Both sets have their 
virtues but they cannot be identical due to the state-dependence of the SIC 
mean field. This suggests to introduce throughout a second basis set ifi which 
is related to the set ip a by a unitary transformation within occupied space, i.e. 

N 



lpa = ^2 i Pi U ia ■ (4a) 



8=1 



We exploit the freedom of a unitary transformation Ui a to choose a basis 
which diagonalizes the constraint matrix \p a . This then yields the diagonal 
SIC equations 

hsicM = e%\<Pi) ■ (4b) 

together with the symmetry condition 

u ia : = {ipp\Up - U a \tp a ) , (4c) 

which now serves to determine the coefficients Ui a of the transformation (j4al) 
for given ip^. Note that the symmetry condition (j4c|) looks similar to Eq. (l3c|) . 
but is used somewhat differently here, namely to define the coefficients of 
the unitary transformation (j4ap . The /isic is defined in detail by the SIC 
equations Q. The new Eqs. (jl]) combine the solution of the former SIC equa- 
tions, yielding the ip a , with the direct evaluation of the energy- diagonal basis. 
What expense is concerned, there is no advantage as compared to the previous 
scheme, i.e. solving first the SIC equations and diagonalizing then the con- 
straint matrix X^ in a second step. But it serves from a formal point of view 
as a preparatory step for SIC-OEP. Although Eq. (l4b| is an eigenvalue equa- 
tion, the corresponding Hamiltonian hsic is obviously non-local, see Eq. ([3d]) , 
which complicates the numerical handling. This is why approximations have 
been searched for since long in order simplify calculations. 



3 SIC-OEP 



3.1 SIC-OEP with double set 



In section 12.21 we have shown that the definition of single-particle energies 
in connection with SIC naturally leads to the introduction of a double set of 
single-particle states connected by a unitary transformation. We transfer now 
this generalization to OEP. The idea beyond SIC-OEP is the following. We 
start from a set if>i, being solution of the OEP equations : 

[k^(T)-V (r)]<p i (r)=e i <p i (T) , (5) 

where Vo is a local and state-independent potential which needs to be opti- 
mized to minimize the SIC energy, Eq. ([1]). Formally, this amounts to param- 
eterize the SIC wavefunctions through Vo, to consider the SIC energy as a 
functional of Vo, and to perform variation with respect to Vq. The \(pi) can 
thus be written as |y?j°) and these are associated to the corresponding set 
of iV'a ) connected by the unitary transformation Eq. (T4"a| . For the sake of 
simplicity, we keep in the following the notation \<pi) for the |y?i°) and use 
correspondingly \ip a ) for \ip\ 



Z°) 



The optimized effective potential Vo(r) is then found by variation 5Esic/8Vo(t) - 
0. One can show that the thus optimized Vq can be written as a sum of three 
contributions : 



V = V S + V K + V C , (6a) 

vfe=£— «< , (6b) 

i p 

Vk = ^2— ! — {<fii\V -Vi\(f>i) , (6c) 

v Cvwm , (6d) 

z i p 

with the density p = J2i \<fi\ 2 - The term Vs is the Slater contribution, Vk + Vs 
constitutes the KLI approximation, and Vq is the remaining genuine OEP 
term. The quantities v% and Pi entering the various contributions are defined 

as 



Vi = ^u* ia -^U a , (7a) 

ft (r) = /dr'{V (r')-^(r')}^G 4 (r,r') , (7b) 

Inserting ( Hal) and ( T7a|) in the quantities involved in ( l6bll6dl) . one obtains : 



lib I 2 
V B = y Z^-U a , (8a) 

a P 

(<fii\V - Vi\(f>i) = ^2u* a U0(^\V o - U a \ip a ) , (8b) 

P l (r) = E<J^'(V (r')-U a (r')) 

a J 

4> a (r')G t (r,r') , (8c) 



We now exploit the interesting property that the ip a are localized [22] ■ This 
means that, at a given r, one single ip a mostly dominates the other wave- 
functions iftp. This amounts to have almost vanishing (18bl) and (IBcj) . and thus 
allows to neglect in Vo the two terms Vk and Vq defined in Eqs. (|6c|) and (J6dj) 
respectively The final OEP result thus naturally reduces to 



em(r) = [h (LDA) (r)-V (r)}&(r) , (9a) 

Vo*E— iWI&l 2 ] > ( 9b ) 

a P 

= (^ /3 |f/ L DA[|^r]-f / LDA[|^r]|^), (9C) 

ip a = J2& U ia , (9d) 

i 

with the symmetry condition ( |9c|) explicitely rewritten for completeness. The 
mean field equation ( |9a|) generates the set </?, of occupied states, while the 
unitary transformation (19dp serves to accommodate the symmetry condition 
(pel) which, in turn, defines the "localized" set ip a that enters the Slater mean 
field Vq as given in Eq. (l9b|) . Note that this equation has the form of a Slater 
approximation [!Mf2T] but is constructed from the localized ipa and applied 
to the (f>i. We thus call this new scheme "Generalized SIC Slater" approxima- 
tion, which differs from the usual Slater scheme because of the two basis sets 
involved here. 



3.2 Recovering "standard" SIC-OEP 



The SIC-OEP scheme presented in section 13.11 differs from usual SIC-OEP 
formulations in that it employs simultaneously two complementing sets of or- 
bitals ip a and (pi. It is instructive to step back to "standard" SIC-OEP. Indeed 
the recovery procedure sheds some light on the relevance of approximations 
performed in standard SIC-OEP approaches. One recovers the standard for- 
mulation by simply setting the unitary transform Eq. (Hal) to unity, i.e., 



Uia -> 8 ia (10) 

which renders the sets i[) a and <pi identical. The symmetry condition becomes 
obsolete because there is no unitary transformation to be optimized. The cru- 
cial Vi then simply reduces to ^lda[|^j| 2 ]5 see Eq. (J2bl) . and the usual Slater 
and KLI approximations are directly obtained from Eqs. ()6aH6cp . KLI means 
to approximate Vq ~ Vjg + Vk while the usual Slater approximation goes one 
step further to Vb ~ Vs, but then applied in a situation with only one set of 
wavefunctions. 

The usual reasoning to validate the KLI approximation, i.e., neglecting the 
correction term (16dp . relies on the oscillating nature of p»(r) which makes it 
small in the average. Going from KLI to Slater is then usually justified only 
in the case of homogeneous or well localized systems (described by (f>i). Our 
Generalized Slater potential, as built from localized orbitals (described by ip a ), 
thus naturally appears as a good and well justified approximation to full OEP. 
It remains to see how it performs in actual computations. 



4 Test for a ID dimer molecule 



As first applications of the generalized Slater approximation, we consider two 
test cases. First, we consider a dimer molecule in one dimension (ID) which is 
a widely used schematic model for critically probing crucial structural and dy- 
namical features of an approach, see e.g. [2H]. In particular the studies of [2H] 
have shown that ID models carry all features of (de-)localization. The two elec- 
trons in the dimer have the same spin in order to test the exchange term and 
the self-energy error. As an interaction, we use a smoothed Coulomb poten- 
tial, it>ij = e 2 /J(xi — Xj) 2 + a?, where the parameters ay for electron-electron, 
electron-ion and ion-ion interactions are tuned to reproduce typical molecular 
energies. Taking this interaction, we develop within LDA an energy functional 
for the exchange term. Working at the level of exchange only allows to have 
Hartree-Fock (HF) calculations as a benchmark to which approximations can 
be compared. The second test case will be naturally a "localized" system with 
more electrons, namely a C atom with 4 electrons (3 spins up and 1 down), 
computed in realistic three dimensions. 

Fig. [1] shows the result for the dimer case. One finds in the upper panel the 
dimer binding energy -Ebond = -Earner — 2-E ato m as a function of atomic distance 
for the various approximations. The lower panel shows the dipole polarizabil- 
ity old which is a sensitive test for DFT considerations [30f3~T] . The dimer 
binding energy in Slater approximation (dashed-dotted curve) shows an un- 
natural bump with a maximum at distance of 5 a§. This is caused by a trend 
to derealization in this approximation. The same feature happens in KLI (not 
shown here). SIC instead, which has a tendency to localization, produces a 
reasonable energy curve (dotted line under the full curve) close to the exact 
HF case (full line). It is gratifying to see that the generalized Slater approxi- 
mation (full curve labeled "GSlat" , almost identical to the SIC one) performs 
equally well. The simple-most LDA does also perform fairly well what the 
dimer binding curve is concerned (dashed line almost identical to the SIC 
one). A different perspective is provided by the results for the polarizability. 
LDA is notoriously off the benchmark HF result. SIC performs very well. SIC- 
Slater, however, fails again. But generalized Slater (full curve) does follow the 
SIC result nicely for physically relevant bond distances. 

Fig. [2] analyzes the findings by comparing the terms neglected when stepping 
from OEP down to (generalized) Slater. Mainly because of the oscillating 
behavior of Pi defined in Eq. (I7bl) . the mean values of V« and Vc are very small. 



Thus we plot in Fig. [2] their root-mean-square expectation value, i.e. J{V 2 ) 
with (V 2 ) = f dxV?(x), , normalized by that of Vs. The r.m.s. expectation 
value of the KLI term Vk is generally smaller than that of Vc- The gain in 
smallness by generalized Slater is dramatic (mind the percentage scale). Even 
there, one observes some growth of the r.m.s. expectation value of Vc for small 
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Fig. 1. Upper : Born-Oppenheimer surface for a ID dimer molecule with two elec- 
trons. Compared are results from HF, LDA, SIC, usual Slater approximation, and 
the generalized Slater approximation (GSlat). Lower : Polarization as a function of 
inter-atomic distance for the five approximations considered. 

distances. This coincides with the deviation in the polarizability observed in 
Fig. CD for very small distances. 



As a test case for a compact electron distribution, we have considered the C 
atom using fully three dimensional calculations [32]. The starting point is a 
calculation with exact exchange and the approximations are compared con- 
sistently with LDA for exchange only. Fig. [3] shows the relative deviations for 
energy and polarizabilities a (note that the electron cloud is slightly asymmet- 
ric such that a depends on the spatial direction). The conclusion is obvious. 
SIC is a good approximation to the exact calculation and generalized Slater 
does equally well, while simple Slater or LDA show larger deviations. 
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Fig. 2. R.m.s. expectation value of the neglected KLI-OEP terms, Vk and Vc, 
rescaled by that of the Slater potential Vjg : in the usual Slater approximation 
(thin curves), see Eqs. (|6cl6d|) . and in our generalized Slater one (thick curves), 
calculated as in Eqs. (|6c|6d|) but with the ij) a instead of the fy. Results for the ID 
dimer molecule with two electrons are drawn versus inter-atomic distance and are 
given in percentages. 
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Fig. 3. Comparison of various approaches for the C atom. Left panel: Deviation of 
binding energy from the Hartree-Fock benchmark value. Right panel: Deviation of 
dipole polarizabilities (in x, y, and z direction) from the Hartree-Fock values. 

5 Conclusion and outlook 



We have presented a formulation of SIC-OEP which employs two different sets 
of N single-particle wavefunctions. One set is taken for the solution of the OEP 
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equations, thus diagonal in energy and most likely delocalized. The other set is 
used in setting up the SIC energy which becomes minimal for localized wave- 
functions. Both sets are connected by a unitary transformation which leaves 
key features as, e.g, the total density invariant. Using that double set allows 
to accommodate two conflicting demands, energy diagonality versus locality. 
The unitary transformation is determined by minimization of the SIC energy 
which leads to what we called the symmetry condition, a key building block 
of the SIC equations. The localized character of the SIC-optimizing set is well 
suited to justify the steps from OEP to KLI and further to Slater approxima- 
tion. We call that scheme generalized Slater approximation. By virtue of the 
double-set technique, it has a wider range of applicability than straightforward 
KLI or Slater approximation. We have tested the scheme on a ID model of 
a dimer at various distances and on the C atom in full 3D, comparing LDA, 
standard Slater approximation, generalized Slater approximation, and SIC 
with full Hartree-Fock as a benchmark. The numerical tests demonstrate that 
the generalized Slater approximation performs remarkably well. The present 
development is to be considered as a first step. A generalization to a time- 
dependent version is obvious, but raising several new consistency conditions 
(energy conservation, stability). Work in that direction is in progress. 

This work was supported, by Agence Nationale de la Recherche (ANR-06- 
BLAN-0319-02), the Deutsche Forschungsgemeinschaft (RE 322/10-1), and 
the Humboldt foundation. 
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